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ABSTRACT 

Hierarchical structure formation implies that the number of subhalos within a dark matter halo 
depends not only on halo mass, but also on the formation history of the halo. This dependence on 
the formation history, which is highly correlated with halo concentration, can account for the super- 
Poissonian scatter in subhalo occupation at a fixed halo mass that has been previously measured 
in simulations. Here we propose a model to predict the subhalo abundance function for individual 
host halos, that incorporates both halo mass and concentration. We combine results of cosmological 
simulations with a new suite of zoom-in simulations of Milky Way-mass halos to calibrate our model. 

We show the model can successfully reproduce the mean and the scatter of subhalo occupation in these 
simulations. The implications of this correlation between subhalo abundance and halo concentration 
are further investigated. We also discuss cases in which inferences about halo properties can be affected 
if this correlation between subhalo abundance and halo concentration is ignored; in these cases our 
model would give a more accurate inference. We propose that with future deep surveys, satellite 
occupation in the low-mass regime can be used to verify the existence of halo assembly bias. 

Keywords: dark matter — galaxies: halos — methods: analytical — methods: numerical 


1. INTRODUCTION 

Bridging our understanding of the processes of galaxy 
formation and of the evolution of dark matter halos re¬ 
mains one of the primary challenges in modern cosmol¬ 
ogy. While 7V-body simulations provide detail about the 
formation and evolution of dark matter halos, it is still 
observationally challenging to directly probe their prop¬ 
erties. Nevertheless, extensive work over the past decade 
has used observations of galaxy’s spatial distributions to 
constrain models of the galaxy-halo connection, which 
reveals how galaxies form in halos (e.g., Berlind & Wein¬ 
berg 2002; Zehavi et al. 2011; Reddick et al. 2013). As 
new observations become more precise, it is crucial to 
understand possible systematic uncertainty and bias in 
those models. 

The two main characteristics of a dark matter halo are 
its mass, usually calculated by setting a spherical over¬ 
density region, and its formation history. The latter is 
also highly correlated with the density profile of the halo, 
and hence with the concentration and with the maximal 
circular velocity 'Umax of the halo (Wechsler et al. 2002). 
Halos of the same mass but different formation history 
can have very different characteristics or reside in differ¬ 
ent environments (e.g., Bullock et al. 2001; Allgood et al. 
2006; Maccio et al. 2007). 

The abundance of subhalos within a dark matter halo 
most strongly correlates with the mass of the halo (e.g., 
Kravtsov et al. 2004). Nevertheless, at a fixed halo mass, 
the subhalo abundance also correlates with the formation 
history of the halo (Zentner et al. 2005; Zhu et al. 2006; 
Ishiyama et al. 2009). This correlation, despite its signifi¬ 
cance in modeling satellite occupation, is often neglected, 
mostly because it does not manifest itself when the Pois¬ 
son scatter is comparable to the number of subhalos in 
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consideration. Satellite occupation, or richness, is often 
used as a proxy of host halo mass, especially for optical 
observations of clusters (Rozo et al. 2009, 2010). The 
scatter in the mass distribution inferred from richness 
can be underestimated if this correlation with concen¬ 
tration is neglected. 

In this work, we investigate again the correlation be¬ 
tween subhalo abundance and halo concentration, and 
propose a simple model that describes this correlation. 
This model can also be used to extend the subhalo abun¬ 
dance function for a given host halo beyond the resolu¬ 
tion limit, and enables us to evaluate how this correlation 
may manifest in a range of observable statistics. 

The simplest approach to extend the subhalo abun¬ 
dance function beyond the resolution limit is to extrap¬ 
olate a parametrized subhalo abundance function. The 
subhalo abundance function is most commonly modeled 
by a power law, and the parameters of the model can be 
calibrated against simulations. Studies have shown this 
method describes the subhalo abundance functions in N- 
body simulations very well (Gao et al. 2004; Kravtsov 
et al. 2004; Giocoli et al. 2008; Springel et al. 2008; An¬ 
gulo et al. 2009; Boylan-Kolchin et al. 2010; Ishiyama 
et al. 2013; Cautun et al. 2014b), at least for host halos 
in a narrow mass range. 

In order to calibrate this kind of model over a wide 
range of mass, usually a suite of cosmological simulations 
and zoom-in simulations is needed. Zoom-in simulations, 
such as the Aquarius and Phoenix simulations (Springel 
et al. 2008; Gao et al. 2012), are particularly powerful for 
measuring subhalo abundance function at high resolution 
but still with reasonable costs. However, if one wants to 
study the halo-to-halo scatter in the subhalo abundance 
function, a fairly large sample size is required. More 
recently, two re-simulation suites have been completed 
with tens to hundreds of simulations in specific small 
mass ranges: the Rhapsody (cluster-mass halos, Wu et al. 
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2013) and ELVIS simulations (Milky Way-mass halos, 
Garrison-Kimmel et al. 2014). 

While these fitting models can usually describe simula¬ 
tions fairly well, they often capture the minimal relevant 
physics for the particular questions that are being ad¬ 
dressed. A more elaborate approach is to consider the 
assembly histories of dark matter halos and the evolu¬ 
tion of halo mass function (Yang et al. 2011). One can 
further consider more relevant subhalo dynamics when 
modeling subhalo abundance beyond the resolution limit 
by tracking the orbits of subhalos and adding subhalos 
that do not appear or are disrupted in simulations (Zent- 
ner et al. 2005; Jiang & van den Bosch 2014; van den 
Bosch & Jiang 2014). Instead of fitting the abundance 
function, this kind of approach considers most physical 
details, but at the same time can be more difficult to 
constrain. 

In this work, we focus on an empirical model which 
directly uses mass and f max of the host halo to pre¬ 
dict the subhalo abundance function, and calibrate the 
model against cosmological and zoom-in simulations. 
This model is essentially the simplest possible model of 
subhalo abundance function that takes halo formation 
history into account. In principle, a more sophisticated 
model (i.e., models that track subhalo evolution) could 
produce similar results. However, our simple model pro¬ 
vides a straightforward way to evaluate this correlation 
between subhalo abundance and halo formation history, 
and to evaluate its implications for various observables. 

This paper is organized as follows. In Section 2 we 
describe the simulations used in this study. In Section 3 
we first discuss the correlation between subhalo abun¬ 
dance and halo formation history, and then we describe 
and calibrate the model which predicts the subhalo abun¬ 
dance. In Section 4, we further discuss the implications 
of this correlation between subhalo abundance and halo 
concentration. We summarize this paper in Section 5. 

2. SIMULATIONS 

In this study we use a cosmological simulation 
c 125-2048 and also present a new set of zoom-in sim¬ 
ulations of Milky Way-mass halos. 

The c 125-2048 box 3 is a dark matter-only cosmologi¬ 
cal simulation run with L- Gadget (based on Gadget- 
2, Springel et al. 2001; Springel 2005). The box has 2048 3 
particles and a side length of 125 Mpc h -1 , resulting in a 
particle mass of 1.8 x lO 7 M 0 h _1 . The softening length 
used is 0.5 kpch -1 , constant in comoving length. The 
cosmological parameters are = 0.286, = 0.714, 

h = 0.7, erg = 0.82, and n s m 0.96. The initial condi¬ 
tions are generated by 2LPTIC 4 (Crocce et al. 2006) at 
z = 199, with the power spectrum generated by Camb. 5 

The new suite of zoom-in simulations consists of 46 
Milky Way-mass halos, selected from the c 125-1024 box 
(see footnote 3), which is a low-resolution version of the 
c 125-2048 box. The parameters and initial conditions 
of these two boxes are identical, but c 125-1024 contains 
only 1024 3 particles and starts at z = 99. All the selected 
halos fall in the mass range M vir = 10 121±0-03 M o in the 

3 Provided by Matthew Becker (M. Becker et al. 2015, in prepa¬ 
ration) 
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c 125-1024 box. The initial conditions of these zoom- 
in simulations are generated with the publicly available 
Music code 6 (Hahn & Abel 2011), and are matched to 
the cosmological box up to the 1024 3 scale. The La- 
grangian volume where the highest-resolution particles 
are placed is set by the rectangular volume which the 
particles within 10R v ir of the present-day halo occupied 
at z = 99. The mass of the highest-resolution parti¬ 
cles in the zoom-in simulations is 3.0 x lO 5 M 0 h _1 . The 
softening length in the highest-resolution region is 170 
pc h -1 comoving. Figure 1 shows the images of 6 of these 
zoom-in simulations. Figure 2 compares the concentra¬ 
tion distribution of this sample of Milky Way-like halos 
with the full sample in the mass range in the c 125-2048 
box. The concentration distribution of the selected sam¬ 
ple is slightly wider than that of all the host halos in the 
mass range. 

In the analysis, we use Rockstar 7 for halo finding and 
Consistent Trees 8 for tree building (Behroozi et al. 
2013a, b). The halos are defined with A vir ~ 99.2 for 
this cosmology. Subhalos are defined as halos that are 
within R v i r of any other larger halo. Halos that are not 
a subhalo are called host halos throughout this paper. 

The particle mass of a simulation cannot be directly 
translated into the maximal circular velocity, tWx, to 
which the simulation converges. By inspecting the ve¬ 
locity function, we estimate that a conservative lower 
limit for the convergence of the c 125-2048 box is 40 
kms -1 , and that of the zoom-in Milky Way simulations 
is 9 km s -1 . 

3. MODELING SUBHALO ABUNDANCE 

In this section, we present a framework to model the 
subhalo abundance of individual host halos. We first 
discuss the correlation between subhalo abundance and 
host halo concentration, and observe qualitatively how 
host halo concentration affects subhalo abundance func¬ 
tion. We further argue that for a given host halo, the 
number of subhalos is consistent with a Poisson distri¬ 
bution. Then we describe both the framework and the 
specific parameterization of our model, and calibrate the 
model against the aforementioned simulations. Finally 
we briefly discuss the universality of the subhalo abun¬ 
dance function. 

3.1. Dependence of Subhalo Abundance on Halo 
Concentration 

7V-body simulations have shown that the subhalo 
abundance function averaged over a sample of host halos 
of a similar mass approximately follows a power law, and 
its form is nearly universal for different host halo masses 
when scaled properly (e.g., Gao et al. 2004; Kravtsov 
et al. 2004; Boylan-Kolchin et al. 2010). Hence, the sim¬ 
plest model of subhalo abundance is to describe the mean 
number of subhalos, (lV su b), as a function of host halo 
mass only. Although this simple kind of model can pre¬ 
dict the mean number of subhalos at a given host halo 
mass in simulations fairly well, it cannot explain the de¬ 
pendence of subhalo abundance on host halo concentra¬ 
tion, as shown in Zentner et al. (2005). 

6 https://bitbucket.org/ohahn/music 

7 https://bitbucket.org/gfcstanford/rockstar 

8 https://bitbucket.org/pbehroozi/consistent-trees 
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Figure 1. Images of the zoom-in simulations of six Milky Way-mass halos, from our suite of 46 halos. The concentration of these selected 
halos decreases from left to right. 



Figure 2. The cumulative distribution of concentration (in log 
scale) for the zoom-in Milky Way halos (red) and all the halos in 
the same mass range in thee 125-1024 box (blue). 


To see how host halo concentration affects the number 
of subhalos, in Figure 3 we plot the mean number of halos 
(including hosts and subhalos) whose T max (or T pea k) is 
larger than 60 kms -1 (or 75 kms -1 ) as a function of 
host halo mass. We plot this relation for all the host 
halos and for only halos with the highest and the lowest 
25% of concentration in each mass bin. We can clearly 
see that halos of high concentration tend to have fewer 
subhalos, and also see that this is not a small effect, 
especially when the halo halo mass is about 10 12 M©/i - 1 . 
We note that at higher host halo mass, this difference 
becomes smaller because high-mass halos have a smaller 
spread in concentrations than low-mass halos. 

We now take a closer look at how concentration affects 
the subhalo abundance on a halo-by-halo basis for host 
halos of the same mass. In Figure 4, we plot the subhalo 
Vmax function for all the zoom-in simulated Milky Way- 
mass halos. The subhalo T max functions in Figure 4 are 
colored according to the concentration of their respective 
host halos. We observe two prominent features: 

1. All these halos fall in a very narrow mass bin 
(smaller than 0.08 dex), yet there is a significant 
halo-to-halo scatter in their subhalo r ma x functions. 
The halo-to-halo scatter seems to affect mostly the 
normalization of the subhalo T max function, and 
the trend roughly follows the concentration trend, 
which is indicated in colors — darker lines sit lower. 

2. On the log-log plot, subhalo T max functions are 
mostly parallel to one another, especially in the 
regime where 7V su b > 10. This suggests the power- 
law index is roughly a constant from halo to halo. 
Also, for each individual halo, the deviation of the 
abundance function from a simple power law is 
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Figure 3. Number of galaxies , i.e., halos (including both hosts 
and subhalos) with a cut in u max (upper) or in u pea k (lower), as 
a function of host halo mass. The black solid line shows all host 
halos, while the blue dashed line and the red dash-dot line show the 
host halos with the lowest and the highest 25% of concentration, 
respectively. 


much smaller than the halo-to-halo scatter when 
Asub is large. 

In Wu et al. (2013), the authors also find that the num¬ 
bers of subhalos in different T max bins are correlated, es¬ 
pecially when A su b is large. This agrees with our findings 
here. 

This correlation between the subhalo number and host 
halo concentration has been found and discussed in, for 
example, Zentner et al. (2005), Watson et al. (2011). 
This correlation can be understood by the hierarchical 
formation of halos: conditioned on a fixed halo mass, 
halos with higher concentration form early, and subhalos 
in these halos are stripped longer to a lower mass and 
Umax) and some could already be completely disrupted 
and merged with the host. Both effects would result in 
a smaller number of subhalos at a fixed velocity cut. 

3.2. Small-scale Poisson Scatter 

It is also known and shown explicitly by Boylan- 
Kolchin et al. (2010) that the scatter in the number 
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Figure 4. The subhalo u ma x function for the 46 zoom-in sim¬ 
ulations of the Milky Way halos. Each line represents one host 
halo and is colored according to the ratio Vhiax/Wr °f the host 
halo. Darker color reprensents halos of higher concentration (larger 
Vmax/Vvir). The gray band on the left shows the regime affected by 
resolution, where the abundance function bends due to unresolved 
subhalos. 

of subhalos is super-Poissonian when the mean number 
is much larger than 1, The authors argue this super- 
Poissonian scatter is a sum of a Poisson scatter and an in¬ 
trinsic scatter (see also related discussion in Busha et al. 
2011b). 

Here we further claim that the Poisson scatter should 
exist on a single-halo basis. That is, given a host halo and 
its environment, the small-scale variation would result 
in a Poisson scatter in its subhalo abundance. On the 
other hand, the intrinsic scatter (or more precisely called 
the halo-to-halo scatter) is then in principle all possible 
scatter among host halos. 

To verify that the subhalo abundance function is al¬ 
ways subject to this small-scale Poisson scatter when we 
consider a single host halo, i.e., 

(V sub | host) ~ Pois((iV sub | host)), (1) 

we run 13 zoom-in simulations of a single halo, with 
different random seeds for the small-scale modes. All 
these 13 realizations have the same simulation setup 
as described above, and also the same large-scale ini¬ 
tial conditions down to the scale of k ~ 16.4 hMpc -1 , 
which is equivalent to 2048 3 particles in the box. This 
scale roughly corresponds to a host halo mass of 2.5 x 
10 10 M o h _1 , or host H max ~ 50 kms -1 . 

Figure 5 shows cr/crp 0 i s , where a is standard deviation 
and crpois = %J (TV), i.e., the square-rooted ratio of the 
variance to the mean of the number of subhalos, in bins 
of u max of the subhalos. The variance and the mean are 
calculated over the 13 halos of the same large-scale ini¬ 
tial conditions. If the number of subhalos in a given u max 
bin follows a Poisson distribution, this ratio would be 1. 
In Figure 5, one can see that at higher values of u max , 
this ratio is less than 1, which is expected due to the 
constrained large-scale modes. At smaller u max , this ra¬ 
tio approaches 1. Although the sample size is small, the 
typical number of subhalos above u max = 10 kms -1 is 
already more than 200. Hence, if the super-Poissonian 
scatter truly exists at the scales within a single host halo, 
one would expect the ratio to be larger than 1 at small 
%ax, scaling similar to the green dashed line, which in¬ 
cludes the super-Poissonian scatter. This test suggests 
that, for a given host halo (and its environment), the 
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Figure 5. The blue line shows cr/crpois in bins of n m ax, calculated 
over the 13 halos of the same large-scale initial conditions. The red 
bands show the l-cr (dark) and 2-cr (light) confidence interval if N 
follows a Poisson distribution and given that there are 13 samples. 
The green dashed line shows the super-Poissonion scatter (Boylan- 
Kolchin et al. 2010, Figure 8) for comparison. 

scatter of its subhalo abundance is consistent with Pois¬ 
son scatter. The super-Poissonian scatter in a fixed host 
halo mass cannot solely come from small-scale modes, 
and should be a result of the scatter in the host halo 
properties at that fixed mass, combined with dependence 
of the subhalo abundance on these properties. 

3.3. Framework of the Model 

Now we present the framework of our subhalo abun¬ 
dance model. We first outline our model that describes 
the number of subhalo for a given host halo, and the pa¬ 
rameters of the model. Then we further present how 
to relate these model parameters to the properties of 
the host halos. In this fashion, we can clearly separate 
the Poisson scatter in each individual host halo from the 
halo-to-halo scatter. 

Mathematically, we can model the subhalo abundance 
function as a counting process. Here the counting process 
we consider is counting over the proxy variable (i.e., u max 
or M v i r ), not over the physical time. Although the math¬ 
ematical term we used is process , we are not considering 
the physical evolution of the subhalo merging process, 
but only the number of subhalos at a given time. 

Let N(v) denote the number of subhalos whose u max 
(or other proxy, which for simplicity we call v) is greater 
than or equal to v. Note that N(v) is always an integer 
and has the following properties: 

N(v i) > N(v 2 ), if vi < v 2l (2a) 

N(v) = 0, ifu>H cut , (2b) 

where V cut is a scale above which there are no subhalos. 
The value of V cut depends on the host halo. 

We further argue that this counting process is an in¬ 
homogeneous Poisson process. That is, the number of 
subhalos in the interval [vi,v 2 ) follows a Poisson distri¬ 
bution and is independent of the counts in any other 
disjoint intervals. We can write 

[N(v i) - N(v 2 )] ~ Pois(A(i>i, 1 * 2 )), (3) 

and 

fviV ( min (v 2 , V cut ) \" 

\Vo) \ Vo ) ' 


X{v 1 ,v 2 ) = 


(4) 
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where Vo is a positive parameter and n is a negative 
parameter, and both could depend on the host halo. Note 
that the parameters Vo and V cut should have the unit of 
the proxy. For example is the proxy is T max , they should 
have the unit of velocity. If one uses M v j r instead as the 
proxy, they should have the unit of mass. 

The expected number of subhalos whose T max > v is 
then simply 




We note that by introducing the V cu t scale, we do 
not need an additional exponential cutoff in the model. 
The average subhalo abundance function naturally drops 
off at the high end, and resembles a exponential cutoff. 
There are two strengths of this approach. First, the pa¬ 
rameter F cut has a clear physical meaning; no subhalo 
can have T max (or any proxy in use) that is larger than 
V cut - Second, when implementing this model, one does 
not need to worry about the chance of having a subhalo 
with a very large T max . The chance of having such an out¬ 
lier is remote but still finite when using an exponential 
cutoff, while in our model the probability of a subhalo 
with T max > V cu t is zero by construction. 

With our framework, there is a straightforward algo¬ 
rithm to create a set of values which represents the set 
of the subhalo r max values of a particular host halo, 
given a known threshold fthres- This algorithm helps 
to generate a mock catalog of subhalos beyond the res¬ 
olution limit. To generate this set, one first draws one 
random number k from a Poisson distribution of mean 
^(^thres) according to Equation (5), with ^thres being 
the minimal possible T max value in the desired set. Then 
one draws k random numbers X \,..., X & from a uni¬ 
form distribution U( 0,1). The desired set would then be 
{f(Xi),...,f{X k )}, where 


f(x) := K, 


-^(^thres) ’ X ~\~ 


Vut 

Vo 


n- 1/n 


(6) 


is the inverse function of Equation (5). 


3.4. Calibrating the Model 

So far we have introduced three parameters that are 
associated with the host halo: W u t, the largest scale a 
subhalo could have; Vo, the overall normalization of the 
subhalo abundance function; and n, the power-law index 
(log-log slope) of the subhalo abundance function. In 
principle, the values of these three parameters in different 
host halos do not need to follow any universal relation, 
and can depend on any host halo property. Nevertheless, 
since the dark matter halos in dissipationless simulations 
do have many universal properties, it is plausible that 
some universal relations relating these three parameters 
to the host halo properties would already make a good 
approximation. 

For conventional models that describe ( N) as a func¬ 
tion of host halo mass only, one can parameterize the 
variables in Equation (4) as follows 

V 0 = aV viT , (7a) 

V cut =bV vil , (7b) 

(7c) 


Table 1 

Parameter Values 


Proxy 

Redshift 

ao 

a 

fro 

p 

no 

%ax 

0 

0.49 

-0.9 

1.4 

-2.5 

-2.90 

%ax 

1 

0.85 

-1.0 

1.4 

-1.0 

-2.80 

%ax 

3 

1.70 

-1.0 

1.4 

-0.8 

-2.60 

^peak 

0 

0.67 

-0.8 

1.4 

-2.5 

-2.75 


Note. — See Equations (7) and (8) for the definitions 
of these parameters. See text of Section 3.4 for details. 


where Wir refer to the circular velocity at R v [ r of the host 
halo, a, 6, and no are all constants that do not depend 
on any host halo properties. 

However, we already know that the parameterization 
above cannot account for the dependence on halo con¬ 
centration. Here we present a specific parameterization 
that replaces a and b in Equations (7) with functions of 
(Vmax/Wr). Particularly, we set 



(8a) 

(8b) 


where ao, &o> and /? are constant. Here V v [ r and H max 
refer to the host halo, and their ratio can be viewed as a 
proxy of the halo concentration or formation time. When 
a = f3 m 0, this falls back to the conventional model 
which has no concentration dependence. 

With this particular parametrization which incorpo¬ 
rates host halo concentration, we can calibrate the model 
against simulations. With the c 125-2048 box, we find 
the values listed in Table 1 provide decent descriptions 
to both the mean and the scatter of subhalo abundance 
across a wide range of mass. We also find the values for 
two different redshifts (z = 1 and 3) and for using v pea ^ 
as the proxy. Note that if one use T pea k as the proxy in¬ 
stead of T max , the dependence on concentration is slightly 
weaker (see the values of a in Table 1). 

Figure 6 compares simulations with the prediction 
from this model with the parameters listed in Table 1. 
In the simulations, we bin host halos according to their 
mass, in a wide range of masses (1O 12 -1O 14 M 0 /i _ 1 ), and 
measure the mean and variance of number of subhalos 
whose t’max > 50 kms -1 in each bin. For each host halo 
we also predict the number of subhalos with the model, 
and measure the binned mean and variance in the same 
way as with simulations. Then we plot the relative dif¬ 
ference between the model prediction and the simulation 
as a function of host halo mass in Figure 6. The relative 
difference is defined as SX := X mo dei/^sim — 1, where 
X could be the mean (upper panels) or variance (lower 
panels) of number of subhalos in each mass bin. 

As Figure 6 shows, our model can reproduce the mean 
and variance of the number of subhalos in all mass bins 
very well. We also plot the model with no concentration 
dependence (a = /3 = 0) for comparison. While this kind 
of model can reproduce the mean value, it fails to repro¬ 
duce the variance. Especially for the predicted variance, 
our model successfully recovers the scatter in high-mass 
bins, where a model that depends only on mass or the 
Poisson scatter cannot. For halos of the highest and the 


n = n 0 , 
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lowest 25% concentration in each mass bin, our model 
also fits the simulation reasonably well. 

In this work, we do not focus on refining these relations 
to obtain the best mock subhalo abundance function. In 
fact, the essence of this work is to show that with our 
simple model one can already reproduce most important 
features in the subhalo abundance function. There are 
two main reasons for not pursuing the best-fit model here. 

First of all, the parameterization proposed above is 
not unique. For example, one can substitute the ratio 
^max/hvir that appears in Vo with some generic function 
of concentration /(c), or put in a mass/velocity depen¬ 
dency in n. The parameters can also involve other host 
properties, or even be stochastic (i.e., involving random 
variables). Also, while the parameters provide insight 
on the dependence on concentration, they do not bear 
clear physical meaning and the parameterization choice 
is somewhat arbitrary. 

Second, although simulations do provide constraints on 
the model parameters, these parameters are very degen¬ 
erate and the Poisson scatter of individual halos makes 
it very difficult to tightly constrain the best-fit parame¬ 
ters. Multiple sets of values could give equally good fits 
to simulations, and the choice of the objective function 
(statistics to minimize) would also affect the best-fit val¬ 
ues. The reported value in Table 1 are obtained by fitting 
only the mean and scatter of subhalo abundance in the 
full c 128-2048 box in bins of host halo mass (i.e. to min¬ 
imize the two leftmost panels in Figure 6), yet these val¬ 
ues also provide decent fits to the individual abundance 
function as shown in Figure 7. 

As a result, here we do not give meaningful error bars 
on the parameter values, but rather simply demonstrate 
the model’s capability of reproducing the subhalo abun¬ 
dance functions. Until the statistics of high-resolution 
halos improves significantly, we recommend optimizing 
the fit every time for each specific use case. 

3.5. The Power-law Index 

So far we have been fixing the power-law index (log- 
log slope) to be a constant that does not change with 
halo properties when calibrating our model against the 
c 125-2048 box. This assumption is consistent with pre¬ 
vious studies (e.g., Gao et al. 2012). However, due to the 
resolution limit, low-mass host halos in a cosmological 
box do not constrain the index as well as the high-mass 
halos because the number of resolved subhalos in low- 
mass host halos is smaller and subject to larger relative 
Poisson scatter. As a result, the value of no in Table 1 
is mostly set by those high-mass halos in the box. 

To investigate whether the power-law index is indeed 
a constant, we check if the model would work for both 
the zoom-in Milky Way halos and the high-mass halos 
in the box. In Figure 7 we compare the subhalo abun¬ 
dance function in simulations with that predicted by the 
model. We discover that a constant index which can 
fit the subhalo abundance function very well for cluster- 
size halos fails to fit the abundance function for zoom-in 
Milky Way-size halos. The log-log slope of the abun¬ 
dance function is steeper for Milky Way-size halos than 
for cluster-size halos. 

We emphasize again that this mass trend is difficult to 
detect in a cosmological box due to limited dynamical 
range. As shown in the upper right panel of Figure 7, 


at L’max = 50 kms -1 , both the number of subhalos and 
the scatter are still consistent with the prediction from a 
constant slope. 

Recall that the power-law index also changes with red- 
shift, as shown in Table 1: at higher redshift, the log-log 
slope of the abundance function is shallower. The rela¬ 
tion between the power-law index, host halo mass, and 
redshift is also discussed in Zentner et al. (2005), Watson 
et al. (2011). An intriguing question is then whether this 
redshift trend and the aforementioned mass trend in the 
index have the same physical origin. 

Specifically, we find that we can fit the subhalo f max 
functions of the zoom-in Milky Way halos and of the 
cosmological box simultaneously (see the lower panels of 
Figure 7) if we replace the constant index by this relation, 

n = —3.05 v(M, z)~ 0A , (9) 

where 


S c « 1.686 is the critical overdensity, D(z) is the linear 
growth rate, and a(M) is the squared root of the mass 
variance (at z = 0) with a top-hat filter of mass M. 

Figure 8 shows the relation of Equation (9) and com¬ 
pares it with the constant values of no in Table 1. Al¬ 
though this is not a proof of the validity of Equation (9), 
it indeed demonstrates the possibility that the mass and 
redshift trends in the power-law index have the same 
physical origin. To robustly verify this connection be¬ 
tween n and n(M, z) would require several sets of zoom- 
in simulations of halos of different masses, preferably also 
with different cosmologies. This is beyond the scope of 
this work, but worth exploring as simulation suites ex¬ 
pand. 

4. IMPLICATIONS AND DISCUSSION 

So far we have been focusing on subhalo abundance 
function and its dependence on host halo concentration. 
In this section, we discuss its observational implications. 
While we cannot observe dark matter subhalos directly, 
we can certainly count the satellite galaxies that sit in 
those subhalos. Hence, the subhalo occupation can be 
viewed as a proxy of the satellite occupation, subject 
to the effect of baryons on the subhalo abundance func¬ 
tion (e.g., Cui et al. 2012; Vogelsberger et al. 2014). Here 
we ignore baryonic effects and directly translate the sub¬ 
halo occupation above a certain velocity cut to the satel¬ 
lite abundance at a luminosity threshold by specifying a 
galaxy-subhalo connection. 

The simplest relation between subhalos and satellite 
galaxies is a one-to-one relation, 

N sub (> v) = -/V sat (> L(v)), (10) 

where L(v) specifies the correspondence between veloc- 
ity cut and luminosity threshold by matching their abun¬ 
dance functions. This is commonly known as abundance 
matching (e.g., Kravtsov et al. 2004; Vale & Ostriker 
2004), which has been shown to work fairly well for 
predicting measurements such as the correlation func¬ 
tions (e.g., Conroy et al. 2006; Reddick et al. 2013). With 
this abundance matching scheme, the model we intro¬ 
duced in Section 3 directly becomes P(7V sat |M, c), and 
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M vir M vir [M @ h- X ] M vir 

Figure 6. Relative difference between the model prediction and simulation of mean (upper row) and of variance (lower row) of the number 
of subhalos, in bins of host halo mass. The middle and the right columns show the lowest and the highest 25% of concentration, respectively. 
Blue solid line shows the model we present here. The green dashed line is a model that depends only on host halo mass (i.e., a = (3 = 0). 
The red dotted line shows the Poisson scatter given the mean value in each bin. 



[km/s] 


[km/s] 


Figure 7. Subhalo abundance function in simulations (red) and predicted by the model (blue). The shade of colors represents the 
concentration (Umax/Wr) of the halo: the darker the more concentrated. The two columns show two different host halo masses. The upper 
row uses the model with constant index (n = no), while the lower row uses Equation (9). The model with constant index cannot reproduce 
the subhalo abundance function for zoom-in Milky Way halos (upper left panel). The gray band on the left shows the regime affected by 
resolution. 


it implies that satellite occupation depends on both host 
halo mass and concentration. 

A different, but also widely used approach is to use 
Halo Occupation Distribution (HOD). Instead of speci¬ 
fying the galaxy-subhalo connection, standard HOD di¬ 
rectly models the probability distribution of satellite oc¬ 
cupation at a luminosity threshold as a function of host 
halo mass (e.g., Peacock & Smith 2000; Seljak 2000; Scoc- 
cimarro et al. 2001; Berlind & Weinberg 2002; Cooray 
& Sheth 2002). That is, it specifies P(N sat > L|M), 
and this distribution of satellite occupation does not de¬ 
pend on host halo concentration. Nevertheless, one can 
also generalize the HOD to include the concentration de¬ 
pendence and to specify P(AT sat |M, c). Yet most studies 
constraining HOD assume the sole dependence on mass. 

Abundance matching and HOD also differ from each 


other in how the positions of the satellite galaxies are 
assigned. However, in the context of satellite occupa¬ 
tion, the only relevant difference is whether or not the 
satellite occupation depends on host halo concentration 
(at a given host halo mass). It is clear that subhalo oc¬ 
cupation does depend on host halo concentration, but 
the stochastic process of galaxy formation could dimin¬ 
ish this dependence. Nevertheless, it is also possible that 
Equation (10) is only perturbed, and the concentration 
dependence of subhalo abundance still survives and re¬ 
sults in the concentration dependence of satellite abun¬ 
dance. 

In this section, we assume the simple relation of Equa¬ 
tion (10), and investigate the implications of the correla¬ 
tion between concentration and satellite occupation. We 
compare the different inferences between these two mod- 
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Figure 8. The solid lines show Equation (9): the log-log slope 
as a function of mass at three different redshifts: z = 0 (blue), 1 
(green), and 3 (red). The dotted horizontal lines show the values 
of no in Table 1 of corresponding redshifts. 

els (with and without concentration dependence) when 
using satellite occupation as a proxy of halo mass. Then 
we look at the the possible signal of halo assembly bias 
with satellite occupation. 

4.1. Satellite Occupation as a Proxy of Halo Mass 

Satellite occupation, especially in the cluster-mass 
regime, has been used to probe the host halo mass (Old 
et al. 2014, 2015; Oguri & Lin 2015; Rozo et al. 2015). 
Conventionally, this is done within the standard HOD 
framework, which ignores the dependence of satellite oc¬ 
cupation on host halo concentration. Here we would like 
to investigate the effects of ignoring this dependence. We 
consider the two subhalo models, as presented in Fig¬ 
ure 6: one only depends on halo mass like the standard 
HOD, and the other incorporates the dependence on con¬ 
centration as introduced in Section 3. We then take the 
host halos from simulations and populate them with sub¬ 
halos according to these two models. This procedure is 
repeated multiple times to obtain enough statistics and 
to smooth the Poisson noise. 

Figure 9 shows the joint distribution of the host halo 
mass and concentration at a fixed satellite occupation, 
iVsat (^max > 75 km/s) = 100, in the context of cluster- 
size halos. We see significant differences between the 
inferences from the two subhalo models, with or without 
the dependence on concentration. Although the mean 
value of inferred mass does not differ more than 1 cr, 
the inferred distribution of mass is much wider in the 
case with the dependence on concentration, and also in¬ 
cludes many more high-concentration high-mass or low- 
concentration low-mass halos. 

The difference seen in Figure 9 would be especially 
prominent when the number of subhalos in considera¬ 
tion's large compared to the Poisson noise, i.e., 7V sat ^ 
yj Wat • Thus when estimating the mass of galaxy clusters 
with richness or satellite occupation, one should consider 
including halo concentration in the model, especially in 
cases when not only the mean estimator but also the 
resulting inference is relevant. 

To refine the mass estimator for halos of a fixed occupa¬ 
tion, we then need some independent observable to probe 
halo concentration. We discuss three possible choices 
here. 

1. The radial distribution of satellites. If satellites 




Figure 9. The joint and marginal distributions of logarithmic 
concentration (y- axis) and logarithmic mass (x-axis) of all the host 
halos which have exactly 100 subhalos whose v max > 75 kms -1 . 
The upper and lower parts demonstrate the inference from the 
two models: (1) with only mass dependence (upper) and (2) with 
both mass and concentration dependence (lower). Dotted lines in 
the side panels show the same marginal distribution for the other 
model just for convenient comparison by eyes. Both models are the 
same as used in Figure 6. The number in the marginal distribution 
of logarithmic mass shows a value. 

trace the density profile of the host halo, then by 
the radial distribution of satellites could provide in¬ 
dependent information on host halo concentration. 
By comparing the number of satellites in different 
projected radial bins, one may be able to select 
those more concentrated halos in a fixed-richness 
sample. 

2. The luminosity of the central galaxy. For example, 
the abundance matching scheme of Equation (10) 
matches luminosity with v max or ^ P eak instead of 
M v i r , and results in the dependence of luminos¬ 
ity on concentration. Hence a further selection 
on the luminosity of central galaxy may provide 
a tighter mass distribution (see also Reyes et al. 
2008). R. M. Reddick et al. (2015, in preparation) 
also finds a negative correlation between the cen¬ 
tral luminosity and richness at a fixed halo mass, 
which agrees with trends proposed here. 

3. The magnitude gap. In addition to the concen¬ 
tration dependence of luminosity, the magnitude 
gap between the central galaxy and the bright¬ 
est satellite galaxy can further depend on the host 
halo concentration. For instance, as suggested by 
our model, the parameter V cu t itself has a concen¬ 
tration dependence, regardless how luminosity is 
matched to halo properties. It has also been shown 
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Figure 10. Same as Figure 9, but showing the distributions of 
log ApJlk Ub (y-axis) and logarithmic mass (cc-axis). The ma¬ 
genta dashed line in the lowest panel shows the mass distribution 
when selecting only halos whose “gap” is larger than 2.5. 

in simulations that the gap is correlated with the 
formation history of the host halo, and hence with 
concentration (D’Onghia et al. 2005; Zentner et al. 
2005; Dariush et al. 2010; Deason et al. 2013; Wu 
et al. 2013). 

It has been suggested that selecting on magnitude gap 
can refine the mass distribution of a fixed-richness sam¬ 
ple (More 2012; Hearin et al. 2013; Lu et al. 2015). Here 
we revisit this method by considering the correlation be¬ 
tween occupation (richness) and halo concentration. Fig¬ 
ure 10 shows the distributions of magnitude gap and halo 
mass, for a sample of a fixed occupation (100 subhalos 
whose T max > 75 kms -1 , same as in Figure 9), for the 
two subhalo models. Here the magnitude gap is approx¬ 
imated by logTp°^/Tp|^ ub , and can be translated into 
the actual magnitude map by abundance matching. As 
we already learned, the distribution of halo mass is much 
wider (lower panel) than that from the assumption that 
satellite occupation depends on host halo mass only (up¬ 
per panel). Nevertheless, if we apply a further selection 
on the magnitude gap, selecting only halos with larger 
gaps, we can obtain a sample of halos whose mass distri¬ 
bution is much closer to that in the upper panel. 

This may provide a viable method to obtain a sample 
of halos in a narrower halo mass bin, especially in the 
high-mass regime. It has been shown that selecting on 
magnitude gap can indeed narrow the velocity dispersion 
distribution of the sample (Hearin et al. 2013). As for 
halo mass, it remains to be seen how strong these effects 
are in specific observed samples, but we expect that the 
relative impact of the central galaxy luminosity and the 
magnitude gap could be tested in the near future using 


Figure 11. Same as Figure 9, but showing of all the host halos 
which have exactly four subhalos whose u ma x > 30 kms -1 . 

either lensing or X-ray measurements of large samples of 
optically selected clusters with fixed galaxy number. 

4.2. Satellites of Milky Way 

In the context of Milky Way-mass halos, the num¬ 
ber of subhalos in consideration is much smaller, and 
the Poisson noise of individual halos would dominate 
and diminish the difference between these two subhalo 
models. Nevertheless, in Figure 9 we observe a positive 
correlation between the host halo mass and concentra¬ 
tion for this sample of a fixed satellite occupation. This 
positive correlation differs from the commonly known 
concentration-mass relation (e.g., Navarro et al. 1997), 
and can also been seen when the number of subhalos in 
consideration is small. 

Figure 11 shows the joint distribution of the host halo 
mass and concentration at another fixed satellite occu¬ 
pation, TV's at (T max > 30 km/s) = 4. In this case, the 
marginal distributions of mass or of concentration barely 
differ between the two subhalo models. Nevertheless, the 
predicted correlation between mass and concentration is 
fairly different in the two cases. Without the dependence 
on concentration, a sample of a fixed satellite occupation 
basically corresponds to a sample of halos in a mass bin, 
and the correlation between halo concentration and mass 
inherits the usual, negative, concentration-mass relation 
of host halos. On the other hand, with the dependence 
on concentration, the inferred correlation between con¬ 
centration and mass becomes positive. 

This discrepancy again highlights the need to consider 
this dependence of satellite occupation on concentra¬ 
tion when inferring the mass or other properties of the 
Milky Way halo from satellites (e.g., Busha et al. 2011a; 
Rodrfguez-Puebla et al. 2013b, a; Cautun et al. 2014a). 
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If the inference is not derived completely from simula¬ 
tions but with the help of a subhalo model which does 
not account for dependence on concentration, such as 
the conventional HOD, then one might need to consider 
the effect discussed above when interpreting the results, 
particularly the degenerate correlation between concen¬ 
tration and mass. We also note that recent constraint 
on the mass and concentration of the Milky Way from 
dynamical tracers have a negatively correlated degen¬ 
eracy (Wang et al. 2015), while occupation-based con¬ 
straints will have the opposite degeneracy if the concen¬ 
tration dependence is properly accounted for, as demon¬ 
strated here. 

This dependence on concentration also suggests that 
one should take the concentration of the Milky Way halo 
into account when investigating the tension between the 
population of subhalos in 7V-body simulations and that 
of the observed Milky Way satellite galaxies (e.g., Kauff- 
mann et al. 1993; Klypin et al. 1999; Moore et al. 1999; 
Bullock 2010; Boylan-Kolchin et al. 2011; Purcell & Zent- 
ner 2012). While a Milky Way-like halo is conventionally 
defined by selecting on halo mass only, it is clear that the 
concentration of the Milky Way halo could potentially 
change the statistical significance of the aforementioned 
tension. In a follow-up paper, we further investigate 
these implications of this dependence on concentration 
for the Milky Way and its population of satellites (Y.- 
Y. Mao et al. 2015, in preparation). 

4.3. Observing Halo Assembly Bias 

Given that satellite occupation is a direct observable 
that is correlated with halo concentration, it may provide 
a way to observationally detect the halo assembly bias. 
Halo assembly bias has been shown to exist in simula¬ 
tions; particularly it is found that host halos of different 
formation histories or concentrations cluster differently, 

&h(M,c)^6 h (M), (11) 

where b h is the halo bias function (Gao et al. 2005; Wech¬ 
sler et al. 2006; Gao & White 2007). The question we 
want to address here is whether we can measure 

&h(M,7V sat )/5 h (M), (12) 

and if so, whether it implies the existence of halo assem¬ 
bly bias as in Equation (11). 

Instead of calculating the bias function directly, we use 
the mark correlation function (MCF) to probe the bias. 
The MCF is defined as 

MCF(m,r)= V (13) 

( i,j)es r 

where S r = {(i,j) : |x^ — Xj| G [r,r + dr]}, and fh is 
the mean of nii over i. The MCF of a specific mark m 
shows whether the averaged value of this mark for halos 
in pairs is higher or lower than the averaged value of the 
whole sample. For each radial bin S r , we find all pairs 
of halos whose separation falls in that bin and measure 
the mark of those halos. To accommodate the possible 
large range of the mark values, we use the ranks of the 
mark instead of the actual value for m, normalized by 
the total number of different values. If Equation (12) 
holds, we expect either a positive or a negative excess in 
the MCF of Y sat . 



r [Mpc /h] r [Mpc /h] 

Figure 12. The MCFs of concentration (left) and of satellite oc¬ 
cupation (Umax > 60 kms -1 ) (right), for host halos whose mass is 
within 10 11 and 10 11 ' 4 MqH -1 . The red shaded area shows the 
range of MCF consistent with no correlation within 2-cr. 

In Wechsler et al. (2006), the authors found a positive 
excess in the MCF of 7V sat in the regime above M*, but 
were not able to find a similar signal below M*. To in¬ 
terpret these results, recall that for halos below the typ¬ 
ical collapse mass M*, high-concentrated halos are more 
clustered; for halos above M*, high-concentrated halos 
are less clustered. In the regime above M*, halos in pairs 
are on averaged more massive but less concentrated, and 
both characters give a higher 7V sat . As a result, the excess 
in the MCF of 7V sat comes from a mixed effect of both 
mass and concentration, and hence it is easy to detect 
this excess but would be difficult to distinguish whether 
this signal is really coming from halo assembly bias. 

On the other hand, in the regime below M*, the de¬ 
pendence of the clustering strength on halo concentration 
switches sign, but the dependence of 7V sat on concentra¬ 
tion remains the same: host halos that form earlier still 
have fewer subhalos at a fixed mass. As a result, in the 
regime below M*, halos in pairs are on averaged more 
massive and more concentrated, and these two charac¬ 
ters have opposite effects on 7V sat . If a negative excess 
in the MCF of 7V sat is detected, this signal must come 
from the contribution of concentration, or halo assem¬ 
bly bias. However, in Wechsler et al. (2006), there were 
not enough subhalos resolved in the simulation for the 
correlation between subhalo abundance and halo concen¬ 
tration to manifest itself, and hence this signal was not 
detected. 

We first calculate the MCFs of halo concentration and 
of satellite occupation by selecting all resolved subhalos 
whose t> max > 60 kms -1 in our cosmological box, for 
host halos in a mass range, 10 11 -10 11-4 M 0 /i _1 , and plot 
the results in Figure 12. The result we found here is 
consistent with previous studies: significant bias in con¬ 
centration, but not in satellite occupation. This result, 
however, does not directly answer whether or not the 
satellite occupation can probe assembly bias, because the 
variance in 7V sat can be large. As we argued in Section 3, 

(Y sat |M,c) - Pois((Y sat |M,c)). (14) 

For host halos in this mass range, the number of re¬ 
solved subhalos is typically less than 10, even for a high- 
resolution cosmological box (e.g., with a particle mass 
of 10 7 M Q h _1 ). Despite the correlation between subhalo 
abundance and host halo concentration, the scatter in 
subhalo abundance can wash out this correlation, espe¬ 
cially for host halos with few subhalos, and render the 
bias in subhalo occupation unobservable. 








Subhalo Abundance and Halo Concentration 


11 



r [Mpc /h\ r [Mpc /h] 

Figure 13. Same as Figure 12, but shows the MCFs of 
model-predicted satellite occupation down to 60, 50, 40, and 30 
kms -1 . The corresponding number densities are 0.122, 0.216, 0.38, 
1.03 (Mpc /h)~ 3 . 


To verify our conjuncture that Equation (12) would 
hold for low-mass halos if the typical value of 7V sat is 
large (> 10), one would need a cosmological box large 
enough to measure clustering statistics and with a par¬ 
ticle mass of ^ 10 5 M o h _1 , but this kind of simulation 
is still beyond the reach of current computational capa¬ 
bilities. Zoom-in simulations can easily provide a much 
better resolution, but those do not provide large-scale 
statistics. With our model, we can predict the expected 
number of subhalos (satellites) to a lower velocity cut 
(higher number density), while preserving the depen¬ 
dence on host halo mass and concentration. We then 
can quantify at what velocity cut (number density) we 
can start to observe the bias in subhalo occupation in 
low-mass host halos. 

Figure 13 shows the model-predicted MCF of subhalo 
occupation for four different thresholds, in the same mass 
range of the host halos, 10 n -10 11,4 M Q h~ l . The host 
halos are selected from the cosmological box, and for each 
host halo we re-populate its subhalos with our model. At 
Tmax = 60 kms -1 the result can be directly compared 
with the right panel of Figure 12. Since our model by 
construction correlates subhalo abundance and halo con¬ 
centration (Tmax/Kir), the lack of signal in the MCF at 
Tmax = 60 kms -1 results from the Poisson scatter. Mov¬ 
ing the threshold down to u max = 40 kms -1 we start to 
see a clear negative excess in the MCF. As we discussed 
above, this negative excess must originate from the fact 
that paired halos are on averaged more concentrated, and 
hence have fewer subhalos. 

This negative excess in the MCF would manifest in the 
projected correlation function by lowering the one-halo 
term if the low-threshold data is available. With up¬ 
coming deep spectroscopic surveys, such as DESI (Levi 
et al. 2013), data with low thresholds will be accessible 
in the near future. Figure 14 demonstrates the num¬ 
ber of galaxies in a volume-limited sample from two ex¬ 
emplary surveys, assuming the luminosity function re¬ 
ported in Blanton et al. (2003, 2005). Both surveys have 


Lnax [km/s] 

15 20 25 30 40 50 60 80 100 130 



number density [(Mpc /h) 3 ] 

Figure 14. Expected number of galaxies in a volume-limited sam¬ 
ple as a function of number density (and corresponding halo % a x) 
for two example surveys with different sky coverage (given in square 
degrees). 


a detection limit of m r = 19.5, and their sky coverages 
are 290 and 14,000 square degrees, roughly representing 
the GAMA Survey 9 (Driver et al. 2011) and the DESI 
Bright Galaxy Survey, 10 respectively. With the latter 
survey, a volume-limited sample of a few hundred thou¬ 
sand galaxies with m r <19.5 down to the number den¬ 
sity at 0.4 (Mpc /h)~ 3 would be accessible, and this sam¬ 
ple would be sufficient for a precise measurement of the 
projected correlation function. 

We note that although we assume the simple relation 
of Equation (10) in this discussion, this signal has the 
advantage that it is less sensitive to the details of the 
galaxy-halo relation because it only utilizes the number 
of satellites above a certain luminosity threshold, but 
not other properties (e.g., color) of the satellites. Even 
if galaxy formation introduces additional scatter in the 
satellite occupation, as long as this scatter is smaller than 
the halo-to-halo scatter due to halo concentration, this 
signal would survive in the projected correlation func¬ 
tion. 


5. SUMMARY 

In this work, we model the subhalo abundance on the 
basis of individual halos. The framework of our model 
is based on the fact that the scatter in 7V su b for an in¬ 
dividual halo is consistent with Poisson scatter, and the 
additional halo-to-halo scatter in 7V su b for halos in a mass 
bin primarily affects only the overall normalization of the 
subhalo function. For a large sample of halos, we find 
that the subhalo velocity functions of a sample of halos 
in a mass range are nearly parallel to one another. As a 
result, we can model this halo-to-halo scatter by intro¬ 
ducing additional parameters to the model that specify 
the normalization as a function of additional halo prop¬ 
erties. 

We hence present a model which predicts the subhalo 
abundance based on two properties: V v [ r (equivalent to 
mass) and Umax/Kdr (roughly equivalent to concentra¬ 
tion) of the host halos. This model successfully repro¬ 
duces the mean and scatter in the subhalo abundance 
in a given host halo mass bin. It can then be used to 
predict the number of subhalos for thresholds that are 
lower than the resolution limit of the simulation. It also 

9 http://www.gama-survey.org/ 

10 http://desi.lbl.gov/cdr/ 
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enables one to conveniently sample a sequence of f max 
values that represent the subhalos of a given host halo. 

This model further provides plain insight into the de¬ 
pendence of subhalo abundance on halo concentration. 
We found that the halo concentration affects the subhalo 
abundance function mainly through the overall normal¬ 
ization (Vo in our parameterization), but also through 
the “cutoff” scale (V^ ut ). A constant power-law index 
(n) fits the cosmological simulations well; however, we 
also find that an index that depends on halo mass would 
fit the zoom-in Milky Way halos better. This depen¬ 
dence on mass may have the same physical origin as the 
dependence on redshift. 

With this model, we then investigate the observable 
implications of the correlation between the subhalo abun¬ 
dance and halo concentration. We find that when us¬ 
ing subhalo or satellite occupation as a proxy of halo 
mass, one might need to consider using a concentration- 
dependent model, such as the one presented here, to ob¬ 
tain a more accurate inference. We show that ignoring 
this dependence on concentration could result in a biased 
mass inference and an incorrect joint distribution of mass 
and concentration of the sample. Although these biases 
are small, they may become important as other sources 
of systematic errors decrease. 

We further propose that satellite occupation can be 
used to probe halo assembly bias if we can detect all 
satellites which reside in subhalos down to ~ 40 kms -1 . 
Because in the low-mass regime, high-concentration ha¬ 
los are more clustered but have fewer subhalos, this sig¬ 
nal can probe the halo assembly bias in concentration 
and is not degenerate with the contribution from halo 
mass. This method is also less sensitive to the detailed 
galaxy formation processes because it only depends on 
the total count. 
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